if (!require(tidyverse)) { install.packages('tidyverse') }; require(tidyverse)
if (!require(here)) { install.packages('here') }; require(here)
if (!require(raster)) { install.packages('raster') }; require(raster)
if (!require(sf)) { install.packages('sf') }; require(sf)
if (!require(rasterVis)) { install.packages('rasterVis') }; require(rasterVis)
First we need to to load our raster image. This is a multiband raster so we use the stack function from the raster package. For a singleband raster we would just use the raster function.
img <- stack(here('Data', '20210804_162603_13_2419_3B_AnalyticMS_SR_clip.tif'))
Now we can check the raster for some important information.
res(img) # pixel resolution
crs(img) # coordinate reference system (and units of resolution)
names(img) # band names (default)
## [1] 3 3
## CRS arguments:
## +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## [1] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.1"
## [2] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.2"
## [3] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.3"
## [4] "X20210804_162603_13_2419_3B_AnalyticMS_SR_clip.4"
For ease of reference we can give the bands better names at this point.
names(img) <- c('R', 'G', 'B', 'NIR')
names(img)
## [1] "R" "G" "B" "NIR"
Now using plotRGB we can visualize the image.
plotRGB(img, 1, 2, 3)
## Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale): color intensity 324, not in 0:255
Well that didn’t work. The plotRGB function assumes we are using pixel brightness values (16-bit) so it struggles with our surface reflectance values being 256-bit. Here are a few ways to fix that issue.
plotRGB(img, 1, 2, 3, stretch = 'hist') # Histogram stretch
plotRGB(img, 1, 2, 3, stretch = 'lin') # Linear (scaling) stretch
# Both of these will also need stretches applied, but its good to know the scale argument
plotRGB(img, 1, 2, 3, scale = 256^2) # Theoretical maximum value
plotRGB(img, 1, 2, 3, scale = 13108) # True maximum value
Okay, at 3x3 meters, these images are slow to render. To get an extent object for subsetting you could use the drawExtent function. Obviously I can’t do that in a non-interactive session so I’m going to build one from scratch.
aoi <- c(733787.6, 4437162.0, 741327.3, 4443901.0) %>%
matrix(ncol = 2) %>%
extent()
plotRGB(img, 1, 2, 3, stretch = 'lin', ext = aoi)
To speed things up even more we can make a smaller raster by clipping.
img2 <- img %>% raster::crop(aoi)
plotRGB(img2, 1, 2, 3, stretch='lin')
Raster math is pretty straightforward.
ndvi_calc <- . %>% {(.$NIR-.$R)/(.$NIR+.$R)}
ndvi <- img2 %>% ndvi_calc
plot(ndvi)
Now we can plot NDVI as raw values, or set up a logical statement to isolate some values. Now would be a perfect time to use levelplot.
levelplot(ndvi, cuts = 5, margin=list(draw = F))
levelplot(ndvi > 0.85, cuts=1, margin=list(draw = F), colorkey=F)
We can always plot NDVI alone, but sometimes its better to combine the stacks and simply assign it to a color band (green in this case).
img3 <- stack(img2, ndvi)
plotRGB(img3, 1, 5, 3, stretch = 'lin')